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HIGH RESOLUTION FREQUENCY ANALYSIS TECHNIQUES WITH 
APPLICATION TO THE REDSHI FT EXPERIMENT 


I. INTRODUCTION 

The data reduction of the Gravitational Probe Redshift Experiment (GP-A) 
requires the processing and analysis of low-frequency signals with a very hi^ 
degree of accuracy and resolution. This report discusses various investigations 
of computation methods applicable to the data analysis of this experiment. 

The objective of the GP-A is to test the equivalence principle by measur- 
ing the gravitational frequency shift of a hydrogen maser clock in a space probe. 
The frequency of the flight clock is compared with the frequency of an identical 
maser clock in the ground station while the space probe moves up and down 
throu^ the gravitational field of the earth. The frequency comparison is 
accomplished ty transmission of S-Band radio signals [1, 2] . Dr. Vessot of 
the Smithsonian Astrophysical Observatory fSAO) is the Piincipal Investigator. 
The experiment is scheduled for launch from Wallops Island in May 1976. 

The ’’output signal” of the experiment is the beat frequency between the 
ground maser and the flight maser. This signal, which has a frequency of 
approximately 1 Hz, is recorded in digital form during the mission. The goal 
of the experiment is to test the equivalence principle to an accuracy of approxi- 
mately ±50 X 10"®. To achieve this accuracy, the Redshift (gravitational 
frequency shift) contained in the beat signal must be determined with an accuracy 
of approximately ±1 x 10"®, (This figure includes a safety factor of 5 for errors 
in data processing. ) 

The frequency shift as a function of flight time predicted by theory is 
shown in Figure 1, The maximum fractional frequency shift is very small: 

Af/f ~ 5 X 10"*®. However, in terms of resolution required, one has to deal 
with a nonstationary process, and local frequency measurements are meaning- 
less under these conditions. 


The experiment data will be analyzed by comparing flight data with the 
theoretically predicted frequency shift derived from trajectory information 
obtained by tracking. According to theory, the Redshift has a linear dependence 
on the gravity potential difference A <^> (between the clock locations) expressed 
by Af/f = q; A^/c^. The experiment will determine the value of o = 1 ± e to an 
accuracy of ±e to test the linearity of this relationship ( Fig. 2). Dr, Vessot 
has presented an overall data reduction scheme that is essentially based on 
signal processing in the time domain (processing of the phase of the beat 
frequency) « This report investigates several methods of signal processing in 
the frequency domain which could be incorporated in the data reduction scheme 
or could be used as an alternate method to analyze data from the experiment. 
This alternate approach to data reduction would increase confidence in these 
results, in view of the extremely high requirements for experiment accuracy 
and the situation of a one-shot experiment. 

Achieving the accuracy goal of the experiment requires not only careful 
analysis and design of flight hardware systems, but also of data reduction and 
computation methods. The computation methods [e,g, , Fast Fourier Trans- 
form (FFT), Discrete Fourier Transform (DFT), and matched filter! dis- 
cussed in this report are not new. However, applications of these methods with 
the resolution and accuracy requirements of the Redshift Experiment could not 
be found. The main objective was to investigate computation techniques with 
respect to the accuracy, resolution, and computer effort necessary for applica- 
tion to the Redshift Experiment data analysis. The IBM 360/65 computer was 
used for the test runs and simulated signal processing. The key to data reduc- 
tion in the frequency domain is the adaptive filter technique implemented by the 
FFT, Samples of computer printouts have been included to show the small 
numbers and small changes involved in these computations. 


II. GENERATION OF INPUT DATA 


For the development and testing of computation methods , input data 
similar to the expected flight data had to be generated. 

During the ejq>eriment, the beat frequency of the two masers ( -- 1 Hz) 
will be processed by a 14-bit A/D converter with a sampling; rate of 100/s and 
recorded in digital form on magnetic tape. To generate input data t^es, the 
frequency shift expected durii^ flight ( Fig, 1) was computed using the 
analytical expression: 
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where GM/rj ^2 is the gravity potential at the locations of the clocks, Vj and 
V 2 are the velocities of the clocks in an inertial earth centered coordinate frame, 
c is the speed of light, and A includes acceleration terms. A predicted fli^t 
trajectory was used for the computations of tiie frequency shift according to 
equation ( 1) . The first term in equation ( 1) represents the pure Redshift; the 
second term is the second order Doppler effect. (The first order Doppler effect 
is cancelled out by a special scheme of signal transmission and signal processing 
in the ground station. ) 

Three input data tapes covering three segments of the total flight time 
were prepared by Data Systems Laboratory using the UNIVAC 1108: 

biput Data Tape 1 covers flight segment 500 to 3500 s. 

Input Data Tape 2 covers fli^t segment 3500 to 6500 s. 

Input Data Tape 3 covers flight segment 6500 to 9500 s. 

The following information is contained on these tapes: 

a. Flight time in seconds. 

b. Gravitational frequency shift (pure Redshift) in Hertz. 

c. Second order Doppler frequency shift in Hertz. 

d. Additional frequency off-sets and variations to simulate hardware 
characteristics (for later implementation) , 

e. Combined frequency from 2, 3, and 4 in Hertz. 

f. Aspect angle during flight. 

Computations were made in double precision and transposed to tapes which are 
compatible with the IBM 360/65 computer on which all following evaluations 
were carried out. Each word on these tapes can be truncated in computations 
that correspond to the planned 14-bit A/D converter or to other bit numbers in 
the practical range of 8 to 16, 
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The sampling rate of the frequency is 100/s. There are 10 files per 
tape; each file contains 300 s of flight time. There are 10,000 time increments 
per 100 s interval corresponding to an array size of 10,000 in some of the 
( following computations, A total useful experiment flight time of 9, 000 s is 

i equivalent to 900,000 time increments. 


III. ANALYSIS BY PIECEWISE POLYNOMIAL FIT 


For the frequency analysis of the Fledshift Experiment, an avers^ing 
time interval of 100 s is considered optimum on the basis of the noise behavior 
and statistics of the hydrogen maser clock. As the cloclc moves up and down 
through the gravity field, the beat frequency changes continuously. This change 
is rather large compared to the resolution required (10”®) , as shown in Figure 
3. 


A piecewise polynomial curve fit to the frequency shift curve was tested 
as a first approach. This was done by taking the 10,000 data points in frequency 
and calculating: 

a. The least square polynomial fit by orthogonal polynomials 

b. The polynomial coefficients 

c. The fitted values for each time point ( 100 s interval) 

d. The residuals between actual and fitted segments 

e. The mean and the standard deviation of residuals. 

Figures 4 and 5 show typical computer printouts for two different segments of 
flight time for the coefficients of the polynomial giving a beat frequency: 

f = a^ + ajt + 32t^ + ast^ + a 4 t* + ast^ . 

Also, the mean and standard deviation for increasing order of polynomial fit 
are given in Figures 4 and 5, Central Processor Unit (CPU) time for all fits 
is approximately 15 minutes. 
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The result shows that the mean of the residuals decreases with increasing 
order of polynomial fit (as one would expect) , but the spread of the residuals 
increases. It is therefore not possible to match the frequency curve by a poly- 
nomial fit with the required accuracy. The polynomial £q>proach is not converg- 
ing because of the dynamic behavior of the frequency change \shich is described 
by the rela ively complicated analytical function for the frequency shift [equation 
( 1 ) 1 . 


!V. SPECTRAL RESOLUTION AND UNCERTAINTY PRINCIPLE 


The uncertainty relation [3J states that the equivalent widths of a function 
and its transform are reciprocals: 


Equivalent duration x equivalent bandwidth = 1 . 


This implies that if a stationary Redshift signal of 1 Hz had to be observed in 
real time to a resolution of one part in 100,000, the experiment time would have 
to be 10® or 27.7 h. In the digital data pnalysis, a required spectral resolution 
exceeding the limits imposed by the observation time can be obtained by 
artificially extending the original data sequence with zeros. 

To test this technique, a finite sinusoidal pulse corresponding to a fre- 
quency of 1 Hz was anayzed with the FFT and DFT. The required computational 
steps are shown by the how diagram of Figure 6. The sharpening of the spectral 
line, i. e. , the increased resolution of the finite length pulse extended by adding 
zeros, depends on the following parameters according to the relationship between 
the time and frequency domains: 

N — The number of samples taken to represent the time function; this 
is the data block size, normally to base 2 in the FFT [4,5] , 

At — The sample interval, i. e. , the time between samples. 

T — Total length of record T = NAt . 

N/2 — The number of frequency points in the spectral representation; 

this is half the data block size because the frequency information 
is broken into two displays of real ( magnitude) and imaginary 
(phase) part. 
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Af — Frequency resolution or the number of Hertz between frequency points 
( spectral line? ) . 


^n- ax ” maximum frequency of the spectral display. 


Time and frequency domains are related by Shannon’s theorem. 


At 


2 F 


max 


f 


and the frequency resolution is giver, by the total length T of the recorded time 
function, 


Af = 


T 


Therefore, changing one of these parameters will change the other quantities. 

The signal to be analyzed for the Redshlft Experiment is sampled at a 
rate of 100/s by the A/D converter, which means At= 0.01 s. In this case, 
the frequency resolution Af of the spectral analysis for different data block 


sizes N is as follows: 


N = 2‘® = 65,536 

Af = 0.0015 Hz 

N = 221 = 2,097,152 

Af = 0.000048 Hz 

N = 2^2 = 8,388,608 

Af = 0.000012 Hz 

A computer printout of the spectral analysis representation for the case 
N = 2^® and f = 1 Hz is shown in Figure 7. Figure 7 represents only a small 
poilion of the printout which contains 32,768 spectral elements. The frequency 


points ( spectral lines) are given by the letter R displaying the sin x/x shape 
in discrete form. The separation between t^o adjacent spectral lines 
corresponds to 0,0015 Hz. The maximum amplitude (4930) occurs at the 328th 
spectral element. The computations were performed with an IBM 360. The 
maximum length of the array size of 2*® that can be stored in the computer will 
set a limit to the total length of record (l.e, , frequency resolution) in-core 
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computation. With additional disk memory, much larger arrays can be trans- 
formed at the expense of computing time and obtaining superfluous information. 
However, in high resolution spectral analysis, computation time and effort can 
be reduced by computing spectr J elements over a selected limited range of the 
spectrum [61. A flow diagram of this process is shown in Figure 8. A com- 
puter routine scans the ou^ut vector of the FFT to find the spectral element of 
maximum magnitude, A three-point Lagrange interpolation within three spectral 
elements is then applied and a new block size N is selecled to compute fewer 
spectral elements with higher resolution using the D FT [5). 

Typical results of such a selective computation are represented in 
Figure 9 showing the occurrence of the maximum magnitude of the spectral 
elements, A block size of N = was selected for the FFT (corresponding to 
a 100 s beat frequency interval) . The DFT for 50 spectral lines with N = 2*® 
was used to obtain the high resolution of Af = 0,000012 Hz. The fastest and 
most versatile version of the FFT was used on the IBM 360 computer (7J . 

Typical runs on the machine are characterized by in-core mode with 650 K-bj'te 
memory required. CPU time was approximately 2 min, including 7 s for the 
65 K FFT. 

The technique was applied to certain segments of the simulated Redshift 
signal taiJes. The results for two segments with different frequency change rates 
are shown in Figure 10. The 500 to 600 s interval is at the beginning of the 
trajectory and the 3500 to 3600 interval is near apogee. 

The conclusion obtained from these tests is that the FFT and DFT com- 
bined will resolve a stationary signal or pulse train to some parts in 100,000 
when diC record length is extended with a suti'icient number of zeros. The 
methods could be of value to analyze quasi- station ary signals such as those 
caused by the spinning motion and dynamics of the spacecraft. In these cases, 
however, the utmost resolution would not be required. During a short interval 
at apogee, the frequency change will be slow enough to use the FFT and DFT 
methods to determine the beat frequency with high resolution. Except for this 
short period, the Redshift signal cannot be extracted directly by spectral 
analysis bccau^^ ^ the nonlinear frequency change with time prohibits an interpre- 
tation of the spectrum. Consequently, a matched filter technique described in 
the succeeding section w'as pursued as tlie next approach for analysis of the 
beat signal. 


V. .MATCHED FILTER TECHNIQUE 


In many signal processing applications, the use of the matched filter 
technique yields optimum results. This technique has been widely used for 
processing radar signals [8,9], A matched filter has a frequency response 
that is the complex conjugate of the spectinm of the input signal. The output of 
the normal matched filter closely approximates the sin x/x form, which has high 
side lobes, the lai^rest of which is s^me 14. 5 dB below the main lobe. Such 
side lobes can be reduced ty’ application of a spectrum weigliing function. A 
cosine bell was chosen to pro\dde time domain Hanning. Tapering is performed 
in a circular manner so that an even function will remain even and side lobes 
are reduced to approximately 40 dB. 

Figure 11 shows a flow chart for a linear frequency modulated matched 
filter test case. The magnitude of the frequency sweep is comparable to the 
frequency change of tlie Hedshift beat signal. Analytically generated FM signals 
were used to check the progiam. 

The signal A(I,1) - cos (a t + jit^) is compared in a matched filter with 

signals of tlic form B(1, 1) =cos a t + r ../it^) • witli tlie templates 

K 

f, , = 0.99995 

fr 2 = 1.00000 
03 = 1.00005 . 

The dc output (zero element) of tlic matched filter indicates the best fit for 
r> 2 = 1.00000, which corresponds to a match of tlic signal to itself. The matlie- 
matics of the single-tone frequency modulation can be found in Hcference 10. 

An example of computer printout data showing the first 100 elements of 

the matched filter out}^ut for the three different o is given in Figure 12, The 

K 

result verifies the correct programming and soundness of approach. The '"PU 
time for a Mock size of 32,708 elements and including seven FFT operations 
was approximately 2 min. Even though it takes the whole memorj' of the IBM 
300/05 computer with 050 K bytes, the tinv required for the iuialysis of the 
complete Hedshift flight data would be only several hours CPU time. In the 
case of 05 consecutive evaluations of lOO s segments along tlie trajectoiy’, total 
CPU time would be approximately 7-h. 


VI. ADAPTIVE MATCHED FILTER ANALYSIS 


The process of matching two signals with varying frequencies optimally 
in data processing is best automated via the Golden Section Search. In the case 
of the Redshift Experiment, one signal would be the beat frequency measured 
during flight and the second signal would be the computed beat frequency obtained 
from tracking data and theory [ 11 ) . 

The Golden Section Search which is an iterative process to find the 
maximum (or minimum) of unimodal functions is based on the iteration by 
Fibonacci numbers which are defined by: 


Y = Y + Y , with Yn = 0 and Y, = 1 
n +1 n n -1 “ * 


For example, 


Yo = 0 Y 4 = 3 


Yg = 21 


Yi = 1 Yg = 5 


Yj = 34 


Y 2 = 1 Ye = 8 


Yio = 55 


Y3 = 2 


Y 7 = 13 Yii = 89 


The sequence of Y numbers represents the strategy for selecting subdivisions 
along the X axis in the range where the maximum of a fimction is expected. 
Figure 13 shows that the Golden Section Search converges fast to a fraction of 
the original interval of uncertainty. For comparison, a progressive step search 
and a search by interval halving are also shown in Figure 13, [The ratio 
betwv en two consecutive numbers of the Golden Section Search approaches the 
golden ratio: (1 + \T5) /2 . j 


An adaptive matched filter process using the Golden Section Search was 
implemented on the IBM 360/65 computer. The corresj- nding flow chart is 
shown in Figure 14, The input tapes (described in Section II) provided the 
simulated beat frequency (Redshift signal) , For the purpose of this exercise. 
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the input signal is split to provide the two signals to be matched. The Golden 
Section Search method changes one of the signals for each run through the 
matched filter and determines the quality of the template fit. The fit is expressed 
by cv = (lie), with e representing the deviation. Starting with an initial off- 
set of e = i 0.00005, the best fit was obtained for ct - 1.000000 (or e = 
0.000000) with very few iterations as expected, {a = 1.000000 represents a 
match of a signal to itscl^.) Ten min of CPU time was required for a 100 s data 
interval containing 10,000 points. The FFT requires 650 K Ityte memory and 
32,768 element block size. The computer runs performed verified the software 
design. In addition, the following tests were performed; 

a. Noise from a random number generator was superimposed on both 
the amplitude and phase of the beat signal. The results showed, as 
expected, that the matched filter provides optimal noise rejection 
[8,9]. 

b. To demonstrate that the 14-bit A/D converter will be sufficient for 
the ground station, the computed amplitude of the beat frequency 
signal was truncated corresponding to 10, 12, 14, and 16 bits. The 
matched filter output remained maximal at dc for ct = 1. 000000. 

c. It is possible to extend the search method to more than one parameter, 
e.g. , for additional functions of frequency or phase shift. However, 
one quickly arrives at insurmountable difficulties for more than two 
parameters, since both local and global optima are entering the 
picture. Search for additional parameters amounts to cross- 
correlation techniques that are best done in the frequency domain 

[5] . For this reason, the method of extending records with zeros 
was maintained. 


VII. APPLICATION TO REDSHIFTDATA REDUCTION 


Figure 15 outlines one of several possible schemes for application of the 
adaptive matched filter technique to the reduction of Redshift flight data. The 
scheme has not been implemented yet. One of the inputs required for a simula- 
tion run consists of the hardware calibration functions that will be generated 
during the test program of the payload. A brief discussion of the content of 
building blocks in Figure 15 is presented as follows; 


Block 1 is the digitized beat frequency (i.e, , amplitude samples) 
measured during the experiment and recorded on magnetic tape. 

Block 2 represents the FFT applied to the signal from Block 1 before 
input to the matched filter. The DFT and regression analysis (2A and 2B) of 
the flight data provide information concerning frequency shifts and modulation 
effects related to payload spin and spin axis motion if present in the original 
signal [12], 

Block 3 , the matched filter, correlates the measured frequency shift 
with the computed expected frequency shift generated by several other blocks 
discussed hereinafter. 

Block 4. the output of the matched filter, generates the error signal e 
to control the Golden Section Search for an optimum match of the two signals 
via Blocks 5 and 6. The e that is applied to the pure Redshift signal only is 
iterated until the best fit is achieved. The remaining e value is a measure of 
agreement between tiieoretical and measured Redshift. 

Block 7 combines the frequency from several sources to generate the 
expected/corrected signal input for the matched filter. 

Block 8 changes frequency to amplitude, and Block 9 applies the FFT. 

Block 10 provides tiie flight trajectory data for computation of the Redshift 
in Block 11 and computation of the second order Doppler effect in Block 12 . 

Block 13 generates frequency correction to be applied to ihe computed 
frequency shifts in Block 7, These frequency corrections are obtained from 
hardware calibrations, telemetry data indicating the actual flight conditions, 
and trajectory information. Corrections which can be applied include magnetic 
field effects, centrifugal force on the maser cavity, zero-g effect on the maser 
cavity, antenna characteristics, ambient pressure change, and others. The 
resolution of the computational process itself is several orders of magnitude 
above the level set as tlie goal for the experiment (e = ± 0.00001), as has been 
tested and demonstrated in previous sections. 
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VIII. CONCLUSIONS 


The frequency shift experienced by the flight maser is a nonstationary 
event during the entire spacecraft flight within the required spectral resolution 
of 50 parts in a million. A regression analysis of computed frequencies shows 
that the frequency changes are nonlinear within the sample interval of 100 flight s 
and that any least square polynomial fit through the computed Redshift will 
deviate incompatibly within the required resolution. 

A similar result was obtained by a spectral analysis using the FFT at 
medium resolution (0.001) and a DFT at high resolution (0.00001). Again, the 
determination of an instantaneous frequency at the required resolution is incom- 
patible with the nonstationary Redshift event. 

Matching a measured ensemble to a predicted ensemble in the frequency 
domain by a matched filter will provide the means to evaluate the experiment 
data to the requix'ed accuracy. The two obvious constraints to this are: ( 1) the 
expected frequency shift obtained from spacecraft tracking and hardware cali- 
brations have to be known to an accuracy better than the desired experiment 
accuracy, and ( 2) the signal tliat represents the beat frequency between ground 
and flight maser (after first order Doppler elimination) must not be degraded by 
the data recording more than dictated by the desired experiment accuracy. 

An adaptive search for the best fit between the comp’ited and measured 
signal will determine to what accuracy the equivalence principle has been 
verified by die flight experiment in the evaluation of c , a measure of the 
similarity between expected and measured results. 

The matched filter in general provides four desirable advantages for the 
data reduction: 

a. The evaluation is in the frequency domain. The transform necessary 
for this is an averaging process that forms essentially sums of products and 
utilizes more information tlian, for example. Events Per Unit of Time (EPUT) 
methods [ 13J . 

b. The matched filter has optimal noise rejection. The noise level in 
the simulated maser beat signal was set in test runs to several orders above 
the expected signal, and the same results were obtained as with zero noise. 
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c. The matched filter output is invariant to a constant phase shift 
between expected and actual signal, A phase shift that is known to be propor- 
tional to time is detected by the output of the matched filter and can be corrected, 

d. The FFT at the matched filter input will aid the overall signal iden- 
tification and separation as a ’’quick-look” means to detect payload spin decay, 
wobble, or other imexpected malfunctions. 

The computing times in the IBM 360/65 computer proved to be modest and prac- 
tical, It takes approximately 7 min CPU time to determine the best fit between 
a simulated ground signal and an e3q>ected Redshift signal for any 100 s flight 
interval to e ^ 0,00005, 


13 





CHANGE IN A F IN UNITS OF CYCLES/c/100 & 




60C.C0 FKEij = C.33ii6S658E-01 



IG 


Figure 4. Typical computer printout for least square polynomial frequency fit, 

600 to 700 seconds flight time. 
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Figure 5. Typical computer printout for least square polynomial frequency fit, 
3800 to 3900 seconds flight time (near apogee). 
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Figure 6. Uncertainty principle test using finite sinusoidal pulse, FFT 
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Figure 7, FFT specti'al analysis of finite sinusoidal pulse (1,0 Hz) 
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Figure 8. FFT, Lagrange interpolation, 
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Figure 9. Part of DFT printout. 



FREQUENCY 


Figure 10. Typical spectra of simulated beat signals 




A(l. 1)-C0S( + 



B (1,1) = COS {u>t*a^Pth 
ct, = 0.99995 
02 = 1.00000 
tt3 = 1.00005 


A (I, J) X B* (I.J) 


Figure 11. Matched filter -nathod flow chart. 
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Figure 12. Matched filter data printout. 



FRACTION OF ORIGINAL INTERVAL OF UNCERTAINTY REMAINING, 


{ 



Figure 13, Convergence of different search methods. 







Figure 14. Adaptive matched filter flow chart, 











FREQUENCY CORRECTIONS 



15. Data reduction flow chart, 
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